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The //-theorem, originally derived at the level of Boltzmann non-linear kinetic equation for a 
dilute gas undergoing elastic collisions, strongly constrains the velocity distribution of the gas to 
evolve irreversibly towards equilibrium. As such, the theorem could not be generalized to account 
for dissipative systems: the conservative nature of collisions is an essential ingredient in the standard 
derivation. For a dissipative gas of grains, we construct here a simple functional M’ related to the 
original H, that can be qualified as a Lyapunov functional. It is positive, and results backed by 
three independent simulation approaches (a deterministic spectral method, the stochastic Direct 
Simulation Monte Carlo technique, and Molecular Dynamics) indicate that it is also non-increasing. 
Both driven and unforced cases are investigated. 


I. INTRODUCTION 

In 1872, Boltzmann published one of his most important papers [I|. This contribution can arguably be considered 
as the effective birth of kinetic theory, a domain pioneered by Bernoulli, Joule and Maxwell to name but a few, and 
that has turned into an active field of research in mathematics [1], physics Q, and engineering Q. The interest of 
reference [l[ is twofold. First, the time evolution of the velocity distribution function /(v,/) for a dilute gas far from 
equilibrium was derived, under the essential assumption of molecular chaos (Stosszahlansatz), that is, assuming the 
absence of correlations between the pre-collisional velocities of colliding partners @. For the sake of the discussion, 
we omit here the spatial dependence of the velocity distribution. Second and remarkably, Boltzmann also introduced 
a functional H{t) = f /log/dv of the instantaneous and time dependent distribution function /, that can be seen as 
a non-equilibrium entropy: the corremonding iJ-theorem states that this functional has a negative production, 
and vanishes only in equilibrium 0, @] ■ If therefore is a Lyapunov-like functional. This explains in this particular 
context how molecular collisions irreversibly lead to equilibrium. Given that the underlying equations of motion are 
time-reversible (with elastic collisions), this finding ignited a long lasting and heated debate, an account of which is 
not the purpose of the present paper, and can be found in [T^ (see also [ll| and references therein for a more 
technical discussion). In essence, the statistical nature of the //-theorem was not fully recognized in the early days. 
Ironically, when dissipative collisions are considered -as e.g. is the case between macroscopic grains for which the 
equations of motion are irreversible [T^ . [l3-. a well defined time’s arrow is present at the level of interactions, but 
no //-like theorem could be derived so far [IJL [l^. 

The aim of this paper is to present strong hints that a simple generalization of Boltzmann’s original functional 
can be constructed, that exhibits monotonous behavior with time and tends to zero. More specifically, we shall be 
interested in the dynamical behavior of a spatially homogeneous granular gas, made up of a collection of a large 
number of grains undergoing dissipative collisions. The grains are treated as inelastic hard spheres, see e.g [T^.[T^ for 
reviews of the rich phenomenology of this class of systems. We will consider two cases depending on the dynamics of 
the grains between collisions: the free-cooling case on the one hand, in which the grains move freely and the stochastic 
thermostat case on the other hand, in which the grains are driven by some random external force. In the first case, it 
is known from particle simulations that the system reaches an auto-similar regime in which all the time dependence 
in the onc-particle distribution function goes through the instantaneous temperature (defined as the second velocity 
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moment of the distribution) [iMl- In the second case, a stationary state is reached in which the energy lost in 
collisions is compensated by the energy injected by the thermostat 

We first study an Wparticle model where the dynamics of the particles’ velocities are treated as a Markov process. 
A Lyapunov function can be identified exactly implying, under plausible conditions, that the thermostated system 
reaches the stationary state in the long-time limit. In the free-cooling case, the consequence is that all the time 
dependence in the 7V-particle probability distribution is encoded in the instantaneous temperature. This scaling is 
similar to the one proposed in [2lj at the level of the Liouville equation and, of course, implies the one at the level 
of the one-particle distribution. Inspired by the previous model, we propose a Lyapunov functional for the inelastic 
homogeneous non-linear Boltzmann equation that describe the dynamics of the system in the low density limit and 
for N ^ oo [T^. [^. The functional is measured in three independent and complementary numerical techniques, with 
the result that it is always a nonincreasing function of time. These results point in the same direction as those of , 
which focussed on thermostatted systems with emphasis on simplified collision models or kernel. 

The structure of the paper is the following. In section [HI the iV-particle model is introduced and the corresponding 
iL-functional identified. This is then analyzed in section IIIII at the level of the one-particle distribution function to 
in turn define, in section IIVI the candidate Lyapunov function that is a central object in our work. In section IIVI 
simulation results are also shown. Finally, the last section contains some summarizing remarks. 


II. A-PARTICLE DESCRIPTION 


We consider a system of N inelastic particles of mass m, diameter a and spatial dimension d. We will assume 
that a microstate is specified by giving the velocities of the N particles, V(t) = at a given time t. The 

dynamics is Markovian, generated by the following rule: two particles i and j are chosen at random, together with a 
unit, center-to-center vector <t, from which the postcollisinal velocities v', v' follow 


V,- = V,: - 


((T.Vy)cr, 

i + a 


( 1 ) 

( 2 ) 


Here, Vy = — v^-, and a is the coefficient of normal restitution that will be considered velocity-independent. It 

fulfills 0 < a < 1 and collisions are elastic for a = 1. In the hard particle model, the collision frequency is proportional 
to (o’ • Vy). Nevertheless, we will consider a more general model in which this frequency is proportional to (o’ • Vy)'’' 
where 7 is a fixed positive parameter (7 > 0) [ 2 ^. Then, the hard particles model is obtained for 7 = I, while Maxwell 
molecules are for 7 = 0 . Two cases shall be considered, depending on the dynamics of the grains between collisions. 


A. The unforced system 


Let us first address the free-cooling case. With the dynamics specified above, trajectories are generated and the 
state of the system is described in terms of the A-particle distribution function, pyfV, t). The evolution equation for 
this distribution is the generalization for inelastic collisions of the Kac’s equation 


Q ^ 

PAr(V,t) = — ^L(v„Vj)pAr(V,t), 


dt 


i<j 


where we have introduced the operator 

^(v*, Vj)p7v(V) = Vb)'’ 


ai+^ 


ba{i,j) ^ - 1 


PAr(V,t) d(T, 


( 3 ) 


( 4 ) 


and K is a parameter of the model with the only restriction that K (o' • )'>' has dimensions of inverse of time (in 

the hard sphere case it is AT = where n is the density). The operator bs-ii,j)~^ acts on any function of V 

replacing and Vj by the precollisional velocities, i.e. 


b»{i,j) V(V) = /(vi,... ,v,_i,v*,v*+i,... ,vj_i,v*,vj+i,.. .,vn), 


( 5 ) 


found from inverting the law 


1 +a 

V.- = V,,-^-(cr ■ Vy)cr, 


2a 

l-\- a 


CT • Vy)(T. 
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( 6 ) 

( 7 ) 
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It can be seen that Eq. Q admits a special solution in which all the time dependence in the distribution function 
is subsumed in the granular temperature, T{t), defined as 


NT{t)= j dvy ^U,^pAr(V,t). 


By dimensional analysis, this means that 


PN^y jt) = 


VH{t) 


dN 


i=l 


PN 


V 


VH{t) 


( 8 ) 


(9) 


where vnit) = 


2T{t) 


1/2 


is the thermal velocity. In Appendix]^ a consistent equation for (/ijv is obtained and the 
equation for the temperature is analyzed. For 7 > 0, the temperature behaves for long times as 

-2/7 


m « 


for t 


( 10 ) 


where B is a constant. This means that, for 7 > 0, the temperature ‘forgets’ the initial condition in the long time 
limit [ 2 ^ . This important property suggests to work in the following dimensionless variables 


s = K f dt'vQit'), 

Jo 


C = 


V 


where we have introduced an “effective” thermal velocity (proportional at long times to vh) 


vo{t) = 


(i) 


"(a,)- 


( 11 ) 


( 12 ) 


with B an arbitrary constant. Then, the actual thermal velocity is proportional to the effective one in the long time 
limit. The evolution equation for the distribution function in the new variables, 0Ar(C, s), 

cl,N{C,s)=vo{tf^PN(V,t), (13) 

is 

= l^L(c„c,)<^^(C,s) - • CMC,s), (14) 

Kj 


where B = BjK is a dimensionless constant. This equation is equivalent to the one of an inelastic system (with 
a collision rule given by Eqs. ©-©) whose particles are accelerated with a force proportional to the velocity, i.e. 
^ = -fci, and is usually called Gaussian thermostat. 

Let us analyze Eq. (1141) to establish the conditions under which a stationary state is reached in the long time limit. 
General results pertaining to master equations do apply to the probability distribution, (/> 7 v(C,s), see e.g. [l^. Let 
us assume that there exists a stationary solution of Eq. (iMl) . (^|J(C), which fulfills 

l^L(c„c,)</.?‘(C) = |A.c</,?‘(C). (15) 

i<j 


Then, we consider a convex-up function h{x) {h"(x) > 0), bounded from below and defined for a; > 0, from which the 
J^-functional follows as 


M’j^{s) = J dC(j)%{C)h 


4>n{C, s) 

<)>^(C) 


It is shown in Appendix iBl that 

(s) TV — 1 
ds 2 

4'n{C, s) 


J dc J da{ci2 ■ h 


Ih 

UAr(C',s)1 


[ <^?^(C') J 


-h 


h' 


J>n{C',s) 


L ryc) J L J L ryc) <^|*(C0 J 


c^n{C,s) cj}N{C',s) 


(16) 


(17) 
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where we have introduced the notation C' = C 3 ,..., c^v}, and use has been made of the invariance of (/) 7 v(C, s) 

under the change of labels ■<->■ Cj for any i and j. As h(x) is a convex function, the integrand is negative and Jifpf 
decreases monotonically in time. On the other hand, as is bounded from below, it must reach at long times a 
limit in which = 0. In this limit, the distribution is 0Ar(C, 00 ) = lims_j.oo 4>n{C, s) that fulfills 


4>n{C,oo) _ cj)N{C',oo) 


for all C, (T. 


(18) 


Note that, due to the invariance property of 0Ar(C,s) alluded to above, Eq. (fTSl) is also valid when C' is the 
postcollisional velocity vector for any pair of particles, not necessarily the pair 1 and 2. Then, if for any C and U 
such that > t/^ and there exists a sequence of collisions that links C with U, we can 

conclude that (j)N{C,oo) = (()^(C). 

This result is important because a stationary state in the s-variable is related to a scaling of the form given by 
Eq. (|9l) in the original t-variable. Then, if the dynamics of Eq. m is such that, for any initial condition ^^r(C, 0 ), 
the system reaches a stationary state in the long time limit, the dynamics of Eq. m will be such that, for any 
initial condition, pAr(V, 0), the system will reach the auto-similar regime with a scaling of the form given by Eq. (|9|) . 
Moreover, for 7 > 0, this holds independently of the auxiliary parameter B. Eor 7 = 0 the situation is different, but 
it suffices that there exists one B for which the stationary solution of Eq. m exists. 


B. The driven system 


We next treat the case in which, between collisions, the grains are heated by a stochastic force modeled by a white 
noise. This model is referred to as the stochastic thermostat model and has been extensively studied in the literature 
(23 - 1^ . More specifically, the jump moments of the particles’ velocities, due to the thermostat are assumed 

to verify [s^ 


B, 


lim 

At-s-0 At 




eo 


N-1 


l)d,a3 


(19) 


for i,j = 1,.. . ,N and /3 ,7 = 1,..., d. We have introduced the notation Avi^p = Vi^p{t + At) — Vi^p{t), Vi^p{t) being 
the /3-component of the particle i at time t. The parameter is the amplitude of the noise and (... )noise denotes the 
average over different realizations of the noise. The non-diagonal terms are introduced to conserve total momentum. 
Let us remark that, as discussed in Appendix [Cl if total momentum is not conserved, a stationary state is not possible 
(indeed, the center-of-mass velocity follows then a standard Brownian motion, see also HI). The evolution equation 
for the N-particle distribution function, pAr(V,t), is 

^PAr(V,t) = ^^L(v,,Vj)pAr(V,t) -k ^PN^V,t)\no^se, (20) 

i<j 


where the first term of the right hand side of the equation gives the collisional contribution {L(wi, Vj) is given by Eq. 
®), and the second term gives the contribution of the thermostat. If the jumps due to the thermostat are small 
compared with the velocity scale in which pjv(V, t) varies, the usual conditions to derive Eokker-Planck equations are 
fulfilled and the thermostat contribution can be approximated by 


m 


PAf(V,t)|^ 


N d 


d d 


dvi^p dv 


-PN(y,t), 




( 21 ) 


with Bij.p^^ given by Eq. di), independently of the specific probability distribution of the jumps. Then, in this 
limit, the evolution equation for / 57 v(V,t) is 

§iPN(V,t) = ^Y.L(^.,v,)pNiV,t) + eon^)pN(V,t), ( 22 ) 

i<j 


where the operator ^ is defined as 


^(V) 


1 A 1 ^ d d 

2 dv"^ N — 1 ^ dvi 9v,-' 

i=l ‘ i<j 


(23) 
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As in the free-cooling case, it is convenient to work with dimensionless variables. We introduce 

1 

s = Kv:t, C = -, ’ (24) 

the latter having dimensions of a velocity. For the sake of readability, similar names as for unforced systems have 
been employed. In these units, the evolution equation for the distribution function reads 

= s) + e"^(C)<^w(C, s), (25) 

i<j 


pi 

where is the dimensionless amplitude of the noise. 

Performing a similar analysis as in the free-cooling case, and under the same hypothesis, it can be shown that the 
function 




= J dC^%{C)h 


4>n{C, s) 


(26) 


decays monotonically in time and that, for any initial condition, a stationary distribution p^(V) is reached in the 
long time limit. 


III. ONE-PARTICLE DESCRIPTION 


In the previous section, we provided a description at the Wparticle level. Here, we consider the N ^ oo limit case, 
where the problem is expected to admit a closed description in terms of the velocity distribution function, /(c, s), 


(()Ar,i(ci,s) R::!/(ci,s) for A»I, (27) 

that we take normalized to unity, where the j-th marginal of (j)N is defined by 

cj)N,j{ci,... ,Cj) = J dCj+i .. .dcN(l)N{C). (28) 

Integrating Eqs. (iMll and (l25ll over C 2 ,..., cjv and assuming the chaos property 

<(>Ar. 2 (ci,C 2 ,s) RS/(ci,s)/(c 2 ,s), for A>I, (29) 


the homogeneous Boltzmann equation for the two cases is obtained (resp. in the unforced and driven cases) 

^/(ci, s) = J dC2A(ci, C2)/(C1, s)/(c2, s) - y ^ ■ ci/(ci, s), ( 30 ) 

o p ^2 1^2 

—/(ci, s) = J dC2L(ci, C2)/(ci, S)/(C2, s) -f y 'S)- ( 31 ) 

It is worth emphasizing that the above “chaos” notion has been introduced by Kac in order to formalize the idea of 
asymptotic independence of particles in the limit A —>■ oo. 

Until now, no Lyapunov functional for Eqs. O and has been identified. Nevertheless, the analysis made 
in the previous section at the A-particle level suggests the following. Consider a specific example of discussed 
above. Taking h(x) = a; log a;, which is bounded from below by —e~^ for a: > 0, we get 


JfN{s)= J dC</.^(C,s)log 


^Ar(C, s) 




(32) 


In addition, if we assume that the A-particle distribution factorizes for all times in terms of the one-particle probability 
distribution, i.e. 


s) = /(ci, s)... /(cat, s) 

(33) 

0?‘(c) = r‘(ci)...r‘(c^), 

(34) 


and 
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where /®* is the stationary solution of the Boltzmann equation, Eq. (1301) or dSU, the functional becomes extensive 
and transforms into a functional of the one-particle distribution function /(c, s), 


= J dcf{c, s) log 


/(c,s) 


/«*(c)_ ■ 


(35) 


Let us remark that the factorization form given by Eq. (1551) and Eq. (1551) should be understood in the sense of Eq. 
(1291) : it can represent a good approximation at least in the N ^ oo limit. In the elastic limit, the above argument 
can be completely justified and it has been shown that 


lim = / dc/(c,s)log 

N^oo ' 


/(c, s) 




(36) 


a property called as “entropic chaos”. The above limit has been first established for a particular class of well-prepared 


time independent sequences of A^-particle densities (biv by Kac in 
of solutions to the elastic Kac’s equation ([51) in (see also 
(1551) lies in the proof of the convergence 


(see also HI) and more recently for any sequence 
37l|'). The most important difficulty in establishing 


1 

N 


J dCf{ci)...f{cN)\og(l)%{C) ydci/(ci)log/'**(ci) 


(37) 


that one can deduce from a careful use of an accurate version of the central limit theorem. It is worth mentioning 
that the limit (1361) is not the only possible scenario. Consider for instance the A^-particle McKean-Vlasov model 


d 


N 


N 


-<(.jv(C,t) = J2diycMNAC)^N{C,t)) + Y,^cANiC,t) 


(38) 


where the force field term is obtained from the V-body Hamiltonian function Wn and 2-body Hamiltonian 
function a by 


.. N N 

AnAC) = ^o,Wn{C), IVw(C) ^ a(c,-c,)+^|c,|2. (39) 

i,j — l i—1 

It is clear and well-known that the only positive and normalized stationary state is the Gibbs probability measure 
given by 


0^(0 := ^ Zn := / dCe-^-^^\ 

Zn J 


(40) 


The relative entropy J^(s) defined in (1321) is still a Lyapunov functional but now, under some smoothness and 
boundedness assumptions on the 2-body Hamiltonian function a, one can show that the rescaled relative entropy 
J^(s) converges to the free energy, namely 


^J^(s) := J dc/(c,s)log +^Jdcidc2f{ci,s)f{c2,s)a{ci-C2). 


(41) 


This convergence can be rigorously justified mathematically by (1) using the techniques of to prove the propagation 


of chaos on any fc-marginal for this many-particle system, (2) using the technique in 


35 


Section 7] to prove the 


convergence of the rescaled relative entropy. We are interested here in a dilute system - a proviso necessary for the 
validity of the Boltzmann description in which the precise form of the law of interaction between the particles is 
irrelevant, beyond the fact that collisions are dissipative. We do not expect a ‘Hamiltonian fingerprint’ in the limit 
limjv^oo jN^ and we are then led in the next section to conjecture that the relevant form is (1361) . 

Note finally the “gap” between the iV-particle evolution, where the evolution is linear and infinitely many h are 
admissible in the definition (j26L and the mean field limit N = oo where it is crucial to use the extensivity of the 
logarithm function (imposing h{z) = zlnz) so that the relative entropy scales like the number of particles and the 
rescaled relative entropy can converge to an effective relative entropy for the limit nonlinear Boltzmann equation. 
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IV. THE CONJECTURED LYAPUNOV FUNCTION 


We are then in 
relative entropy) [39j 


osition to define our Lyapunov-candidate functional as the Kullback-Leibler distance (also called 
between the time dependent velocity distribution, /(c, s), and its long time limit, /®*(c), 


je{t) = 


dcf{c,t) log 


7(c,0 ' 


(42) 


A convexity argument shows that this quantity is positive [s^, and by construction, it is expected to vanish at 
long times. Our central conjecture is that it does so monotonously in time, i.e. dJif/dt < 0. It is also important 
to emphasize here that with elastic collisions for which the velocity distribution thermalizes and evolves towards 
/®*(c) oc exp(—c^), the above distance J^{t) reduces to the original Boltzmann H{t) functional alluded to above, up 
to an irrelevant constant. It should also be emphasized that Eq. (l42l) is invariant under change of variable c —>■ 
where '0 is some invertible function, an important requirement for an entropy-like functional |42| . 

We first sketch a heuristic argument and then perform numerical simulations. The goal is to prove that the relative 
entropy production is non-negative 


Af) ■= - J dcxdc 2 La{ci,C 2 )f{ci)f{c 2 ) - {I-a) J dc^-^^fic^) > 0 


(43) 


where we keep track of the inelasticity in L^, defined as above. Then, assuming that evolution nonlinear equation 
propagates strong regularity and decay and arguing without full rigor on the functional spaces level, we search for 
minimisers critical points /* of with strong regularity and decay. We then have heuristically 


^a(r)~^i(r) + 0(l-a) (44) 

for a close to 1, and due to entropy - entropy production estimates, we deduce that /* = + 0(1 — a) = -I- 

0(1 — a), see the argument in [l 8 |. Finally by studying the Euler-Lagrange equation satisfied by the minimisers 
f* and performing a perturbative argument in the 0(1 — a) neighborhood, we prove that /* = is the unique 
minimizer locally. The previous argument is perturbative in dissipation. Numerical data suggest, however, that the 
monotonicity of M’ goes beyond the quasi-elastic range. 

We have implemented three complementary and independent simulation techniques to assess and illustrate our 
central statement that dJff /dt < 0: a spectral approach, the Direct Simulation Monte Carlo (DSMC) technique and 
Molecular Dynamics (MD) simulations. We now discuss each method in more detail. 

• In the spectral method the non-linear Boltzmann equation (1301) is directly solved. The velocity distribution 
is truncated, Fourier transformed (assuming periodic boundary conditions), and the evolution of each Fourier 
mode fk is subsequently computed from e.g. a Runge-Kutta scheme. In the driven case, the evolution is given 

by 


M 

Kit)= E -CeMt), 

l,m= — M 


where the so-called kernel modes depend only on 7 and a, and can be precomputed and stored before 
solving the equation, C being a nonnegative constant. This method was first derived in for the elastic case, 
and then extended to the inelastic case in [4l|. It is deterministic and spectraly accurate by nature, preserves 
mass exactly, momentum and temperatures spectraly and costs 0{M^) operations. It is moreover valid for any 
values of a € ( 0 , 1 ]. 

• The DSMC method is widely used in the present context, in aeronautics, and in microfluidics [d^. N particles 
follow a Kac’s walk in velocity space and in the limit of large N^ the corresponding first marginal, /(c, s), evolves 
according to the Boltzmann equation [s^. The method is Monte Carlo in spirit, and thus of stochastic nature. 

• In the MD simulations the exact equations of motion are integrated, starting from a given initial configuration 
of N grains in a finite simulation box of volume, U, with periodic boundary conditions [dj. This method does 
not rely on the putative validity of a kinetic description and by comparing to the outcome of DSMC, provides 
a stringent test of the theory and predictions under scrutiny. In particular, the spatial dependence is fully 
accounted for within MD -unlike in the DSMC approach used where spatial homogeneity is enforced from the 
outset- and does not rely on the molecular chaos assumption. If TV —>■ 00 and in the low-density limit (or more 
precisely, in the Grad’s limit) the first marginal is expected to fulfill the Boltzmann equation. 








Figure 1. Free cooling. Left) (color online) DSMC results for fx with an initial asymmetric velocity distribution made up of 
three sharp peaks. The parameters are a — 0.9 and N = 1000 and the results have been averaged over 10® realizations. The 
distribution is plotted for different values of the number of collisions per particle r. The black solid line corresponds to the 
initial distribution and the bell-shaped red solid line to the distribution at the end of the simulation (r ~ 14). Note that the 
y-scale is logarithmic to better probe the low probability region. Right) Corresponding evolution of The inset shows the 

same data on a linear-log scale. 


In the simulations, the evolution of the one-particle distribution function has been measured for the two models, i.e. 
the Gaussian and stochastic thermostats, using different values of the inelasticity and starting with different initial 
velocity distributions. With that, the functional can be computed through Eq. where the knowledge of the 

late time distribution is required. Hence, cannot be obtained “on the flight”, but is computed after has 
been measured in the simulations. We have taken the grain’s mass, m, as the unit of mass and the initial temperature, 
T(0), as the unit of temperature. In the MD simulations the unit of length is the diameter of the particle, cr. We 
always considered a two-dimensional system of = 1000 disks. The spectral method is used in 2 dimensions of 
the velocity space, with 64 modes in each space directions. It is known that such a number of modes gives a very 
good accuracy, thanks to the spectral convergence of the method. The Gaussian thermostat case has been studied by 
DSMG and MD, while the stochastic thermostat has been addressed via DSMC and spectral methods [45j |. 

Fig. [1] displays DSMC results for a system with dissipation parameter a = 0.90 heated by the Gaussian thermostat 
with B chosen to have unit stationary temperature. The results have been averaged over 10® realizations and the 
initial distribution has been taken asymmetric with three peaks: 


/(v,0) 



ui) + ^( 5 (v - U2) -t- - U3), 

D O 


(45) 


with U 2 


1 

3 


TjO) 

m 


1/2 

(1,1), Ui 


— 3 u 2 , and U 3 = 5 u 2 . In the left side of the figure, the distribution function for Vx, 


defined as fxi'Vxjt) = J dvyf{v,t), has been plotted for different values of the number of collisions per particle, r. 
Clearly, the behavior of on the right hand side is compatible with an asymptotic vanishing for t —>■ 00 , which 

simply indicates that / tends towards /®*. More interestingly, is non increasing, from the shortest times, to the 
largest ones one can reach in the simulations. 

In Fig. [21 a comparison between MD and DSMC results is shown for a system with a = 0.80. The initial distribution 
is the same asymmetric distribution as in the previous case and, again, B has been chosen to have unit stationary 
temperature. The results have been averaged over 5 x 10® realizations in the two kinds of simulations. The density in 
the MD simulations is n = 0.005 (t“^, which corresponds to a rather dilute system. The excellent agreement between 
MD and DSMC is important, not only because it again points to the monotonicity of but also because the 

MD algorithm provides a reference benchmark (“true dynamics”), which does not rely on the hypothesis leading to 
the Boltzmann equation, and in particular does not a priori assume the system to be homogeneous. It should be 
mentioned though that the parameters chosen for MD are such that the system remains in a spatially homogeneous 
state for all times (see e.g. the discussion in Refs [1^ for the free cooling regime). Let us mention that we have 
observed the same qualitative features for a large gamut of initial conditions (symmetric around the velocity origin 
or asymmetric) and different values of the inelasticity in the whole range, 0 < a < 1. 

In Fig. 121 DSMC results are shown for a system heated by the stochastic thermostat. We have considered two values 
of the inelasticity, a = 0.9 and a = 0.95, with an amplitude of the noise, such that the stationary temperature 
is 8.80T(0) for a = 0.9 and 17.13r(0) for a = 0.95. In the two cases, we have started with the same initial flat 
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X 

Figure 2. Free cooling. Evolution of M’it) as a function of the number of collisions per particle for MD and DSMC simulations, 
for an asymmetric initial velocity distribution. Here a = 0.8, N = 1000 and the results have been averaged over 5 x 10^ 
realizations. The inset shows the same data on a linear-log scale, to emphasize the long time behavior. 


H 




Figure 3. Left) DSMC results of the time evolution of for the the stochastic thermostat for a = 0.9, starting with a flat 

distribution. Right) The same but for a = 0.95. 


distribution, in which all the velocities have the same probability in a square centered in the origin in the velocity 
space 


/(v,0) 


a Vx € [-w,w] and Vy G [-w,w] 
0, otherwise 


1/2 , 

. The results have been averaged over 10“ trajectories. Clearly, as in the previous case, the 

functional Jif decays monotonically for all times. Again, as in the Gaussian thermostat case, the same qualitative 
behavior is obtained for other initial conditions and values of the inelasticity. 

Finally, we present in Fig. |4]the evolution of for small normal restitution coefficients, namely a = 0.1 (almost 
sticky particles) and a = 0.25, in the stochastic thermostat case. The spectral scheme is used for these simulations. We 
show our results for both the assymetric distribution (l45l) composed of three peaks (left) and for the flat distribution 
(right). As in the other simulations, we observe in all these cases a monotone decay of the entropy functional 
Thanks to the accuracy of the spectral scheme, and owing to its deterministic nature, we can observe this decay up 
to the machine precision. Although it may be due to numerical artifacts (this behavior can be also observed in the 


with w = 


6T(0) 
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Figure 4. Left) Time evolution of the entropy for the nonlinear Boltzmann equation with stochastic thermostat, solved with 
the spectral scheme (M = 64), for a = 0.25 (solid line) and a = 0.1 (dotted line), starting with the assymetric, three sharp 
peaks. Right) Same, for the flat initial distribution. 


elastic case a = 1), this decay seems to follow two exponential regimes, a very fast one in short time followed by a 
slower one in larger time. Nevertheless, these decays are always exponential. 

All the simulation results point in the same direction: the functional defined by Eq. can be a good Lyapunov 
functional for the free-cooling case (Gaussian thermostat) and for the stochastic thermostat. It is worth emphasizing 
here that considering the naive functional of the type / / log / instead of the Kullback-Leibler distance may lead 
to non-monotonic behaviour, as shown in [^ . 


V. CONCLUSIONS 

In this paper, we have presented strong hints that the functional given by Eq. (1421) can play the role of a Lyapunov 
functional in the context of a dissipative granular gas: it decays monotonically in time, tending to zero in the late time 
non-equilibrium steady state. These results, in agreement with those of Ref [l^, have been shown by three different 
kinds of simulation methods, for a wide class of initial conditions and a wide range of inelasticities, 0 < a < 1. Our 
functional -that takes the form of a relative entropy, or Kullback-Leibler distance- reduces to Boltzmann’s original 
^ in the case of elastic interactions. Its very form can be directly inspired from information theory, where the 
Kullback-Leibler distance plays a prominent role [s^. 

It should be noted that the asymptotic steady state fst enters the definition (H^ . so that we cannot deduce the form 
of the non-equilibrium steady or scaling solution from our functional. This is at variance with the reversible dynamics 
case (elastic collisions, corresponding to a = 1 above), and is the price for losing the energy conservation law. We 
emphasize that we have not been able to prove analytically our central result, which requires to work at Boltzmann 
equation level. The situation is simpler at A^-body level, where a counterpart of the iL-theorem can be shown. We 
nevertheless believe our result is conceptually important; monotonicity of J^{t) with time is a strong statement, the 
derivation of which has been the focus of some effort in the past. 

In the particular free cooling case, we have shown that the A^-particle velocity distribution function reaches an auto¬ 
similar regime in which all its time dependence is encoded in the instantaneous temperature [a similar scaling has 
already been used in the context of the Liouville equation (with spatial dependence) |2l|]. A comparable analysis can 
be put forward for mixtures, binary or polydisperse, where the existence of a scaling solution implies the coincidence 
of the different cooling rates. 


Appendix A: Evolution equation for the temperature in the free-cooling case 

Assuming the scaling form of the main text, Eq. ([^l, the evolution equation for the temperature can be straight¬ 
forwardly obtained by taking the second velocity moment in Eq. (|3]). The result is 


dT{t) 


= -BT{t) 


+1 


dt 


(Al) 




























where we have introduced the time-independent coefficient 
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B = 




4d 


m / 


with 


A = J = J dci y dc2c7/^V52(ci,C2). 


(A2) 


(A3) 


The function ipi is just the integrated scaled distribution, (^ 2 (ci, C 2 ) = / dc^ ... f dc]\[(fi\f{C). Eq. (lAip can be readily 
integrated, leading to the 


m = 


1 , 7-B^' 

_T(0)'y/2 


for 7 > 0 , 


and 


T{t) = r(0)e"-®‘ for 7 = 0 . 

Then, for 7 > 0 the temperature ‘forgets’ the initial condition in the long time limit and behaves as 

-2/7 


m «(4* 


for t 


(A4) 


(AS) 


(A 6 ) 


We point now that the scaling given by Eq. (0 is compatible with Kac’s equation, Eq. By substituting Eq. 
Q into Eq. the following equation for is obtained 


1 dvH{t) d .c^^(c) + l^L(c,,c,)^^(C) = 0, 


Kv]+^ dt dC 


(A7) 


i<j 


where C = NjvHify- As a consequence of Eq. (lAll) . the term (, 4.1 is time-independent, and thus, Eq. (IAtI) 

is consistent. 


Appendix B: Evaluation of 

Before evaluating let us note that Eq. da can be written as a Master equation 

s) = y dUir (C|U)<^^(U, s) - y dVW{V\C)^N{C, s), (Bl) 

with the transition probabilities per unit time, #^(U|C), given by 

ir(U|C) = irc(U|C) + irr(U|C), (B2) 

where 

Wc(UlC) = y X! / - c')(5(uj - c')nfc^ijA(ufc - Cfc), (B3) 

Kj d 

is the contribution due to collision and 

irT(u|c) = -|c. A5 (u-c), (b4) 

the contribution of the thermostat. 

Eollowing the same steps as van Kampen it is obtained that 

d^/v(s) ^ j dC j dVW{C\V)^%{V) {[x{V,s) - x{C,s)]h'[xiC,s)] 

-|-/i[a;(C, s)] —/i[a;(U, s)]} , (B5) 











where 
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c(C,s) = 


^ 7 v(C, s) 


(B 6 ) 


It can be seen that the contribution coming from the thermostat vanishes. Indeed, 

J dC J dUirT(C|U)(^|^(U){[a;(U,s)-a;(C,s)]h'[a;(C,s)]+h[a;(C,s)]-h[a;(U,s)]} 

= 11 dC I dUd(U - C)(^?^(U)U • ^ {[x(U, s) - x(C, s)]h'[x(C, s)] 

+/i[x(C, s)] — h[x{XJ, s)]} = 0. 


(B7) 


The collisional counterpart can be written as 

J dC J dU>rc(C|U)(/.?^(U){[x(U,s)-x(C,s)]/i'[x(C,s)] + /i[x(C,s)]-/i[x(U,s)]} 

= —y— J dC J da- {ci 2 • 0 -)'>'(;i^(C) {/i[x(C',s)] -/i[x(C,s)] +s)][x(C, s) -x(C',s)]}, 

(B 8 ) 

where we have introduced the notation C' = {c'j^, C 2 , C 3 ,..., cat} and we have used that, due to the invariance of 
4>n{C,s) under the change of labels Ci o cj for any i and j, the N{N — l)/2 collisional terms are equal. By 
substituting Eqs (lB7l) and (IB 8 I) into (IB5I) . the expression of the main text is obtained. 


Appendix C: Consistency of Kac’s equation in the stochastic thermostat case 

In this appendix, it is shown that the conservation of the total momentum assumed by the thermostat is essential 
to ensure that Eq. (12^ be compatible with the existence of a stationary state. We assume that a stationary solution 
of Eq. (1^ exists, / 9 ^(V). It then fulfills 

^J2L(^,,v,)p^j^{V)+eoT{V)p%{V)=0. (Cl) 

i<j 

Let us multiply the equation by vf and integrate. The collisional contribution is 

/ dvu?L(v,,Vj)p^(V) 
i<3 d 

= J dSr{or-Vi2y{ba-l){v^ + V2) 

= J dvp%{V) J d&{& ■Vi2y+^, (C2) 


where use has been made of 


(bo- - y{vl +V 2 ) 



Vl 2 )^- 


(C3) 


The thermostat contribution is 


I dvu?r(V)p?*(V) = I J dvp^(V)|^u? = d^o^ 

where the cross-terms do not contribute. The obtained equation is 

A:(iv-1) 


(C4) 


4iV 


dvp?J(V) / dd-(d-• Vi2)'^+^ = dCo- 


(C5) 
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Keeping in mind the above ‘steady-state constraint’, let us now multiply Eq. (ICll) by Vi • V 2 and integrate. The 
collisional contribution is 

/ dvvi ■ V2 T(v„Vj)p?^(V) 
i<j d 

= ^ y dvp^(V) J dSr{Sr ■ vi2)'*'(&ct - l)vi ■ V2 

= ^(l-a^) J dvpJJ(V) j dd-((T • Vl2)'’'+^ (C6) 

where it has been used that 


1 _ q,2 

{b^ - l)vi • V2 = ---(<T • Vi2)^. 


(C7) 


From the conservation of total momentum in a collision, we further have {b^ — l)(vi -|- V2)^ = 0 , so that 


{ba - l)vi • V2 = - 1 )K + '*^2)) 

from Eqs. (IC3I) and (IC7l) . The thermostat contribution is 


(C8) 


Co' I dvvi.V2r(V)p?J(V) 


IV- 1 


. d d 


IV- 1’ 


(C9) 


where only the cross-terms contribute. The obtained equation is equivalent to Eq. (IC5I) but, as indicated above, 
only because of the presence of the cross terms [second term on the right hand side of Eq. (fl^ ]. In other words, it 
is essential that the thermostat conserve total momentum, otherwise, the center-of-mass of the system undergoes a 
Brownian motion in velocity space, incompatible with the existence of a steady state [s^. 
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